2.1 To do

• Check why there are some negative amplicon (probably some matches ahead - Put a condition that amplicon should be >)

2.2 References used

• https://bioconductor.statistik.tu-dortmund.de/packages/3.4/bioc/vignettes/IRanges/inst/doc/IRangesOverview.pdf
• http://kasperdanielhansen.github.io/genbioconductor/

5.1 Compute position on yeast - DO NOT RERUN

# primers <- read_excel(path = "primers_18S.xlsx", sheet ="primers")

# Read from local database
pr2_db <- db_info("pr2_local")
pr2_db_con <- db_connect(pr2_db)
primers <- tbl(pr2_db_con, "pr2_primers") %>% collect()
disconnect <- db_disconnect(pr2_db_con)


ref_seq <- pr2 %>% filter(pr2_accession == "FU970071.1.1799_U") 

# Matches the position of the primers on the yeast sequence
#  - use map2_dfr to get an data frame on output and not a list
#  - use ~ when defining the function

primers_pos <- map2_dfr(primers$sequence,
                    primers$direction, 
                    ~ get_primer_position(.x, ref_seq$sequence, orientation =.y, mismatches = 3))

primers <- bind_cols(primers, primers_pos)

write_tsv(primers, path = "primers_matches.tsv", na = "")

5.2 Build the primer set table 1

# primer_sets <- read_excel(path = "primers_18S.xlsx", sheet ="primer_sets")
# primers <- read_excel(path = "primers_18S.xlsx", sheet ="primers")

# Read from local database
pr2_db <- db_info("pr2_local")
pr2_db_con <- db_connect(pr2_db)
primers <- tbl(pr2_db_con, "pr2_primers") %>% collect()
primer_sets_all <- tbl(pr2_db_con, "pr2_primer_sets") %>% collect()
disconnect <- db_disconnect(pr2_db_con)

if(gene_selected == "18S_rRNA") {
  gene_regions = c("V4", "V9")
  # Just keep the selected primers (V4, V9 etc..)
  primer_sets <- primer_sets_all %>% 
  filter(gene == gene_selected) %>% 
  filter(gene_region %in%  gene_regions) %>% 
  filter(!(primer_set_id %in% 63:66))
  } else{
  primer_sets <- primer_sets %>% 
  filter(specificity ==  "plastid" | (gene == "18S_rRNA" & specificity == "universal" ))
}

primer_sets <- primer_sets %>% 
  left_join(select(primers, 
                   primer_id, 
                   fwd_name=name,
                   fwd_seq=sequence, 
                   fwd_start_yeast= start_yeast, 
                   fwd_end_yeast= end_yeast), 
            by = c("fwd_id" = "primer_id")) %>% 
  left_join(select(primers, 
                 primer_id, 
                 rev_name=name,
                 rev_seq=sequence, 
                 rev_start_yeast= start_yeast, 
                 rev_end_yeast= end_yeast), 
          by = c("rev_id" = "primer_id")) %>% 
  mutate(length_yeast = rev_end_yeast - fwd_start_yeast + 1) %>% 
  select(gene_region, specificity, primer_set_id,primer_set_name, contains("fwd"), contains("rev"),length_yeast, used_in:remark) %>% 
  arrange(gene_region,  fwd_start_yeast)

write_tsv(primer_sets, path = "primers_Table_1.tsv", na = "")



knitr::kable(select(primer_sets,primer_set_id:primer_set_name, fwd_name, rev_name, length_yeast)) %>%
  kableExtra::kable_styling()

6.1 Programing Notes

• Use Biostrings

Accessor methods : In the code snippets below, x is an MIndex object.

• length(x): The number of patterns that matches are stored for.
• names(x): The names of the patterns that matches are stored for.
• startIndex(x): A list containing the starting positions of the matches for each pattern.
• endIndex(x): A list containing the ending positions of the matches for each pattern.
• elementNROWS(x): An integer vector containing the number of matches for each pattern.
• x[[i]]: Extract the matches for the i-th pattern as an IRanges object.
• unlist(x, recursive=TRUE, use.names=TRUE): Return all the matches in a single IRanges object. recursive and use.names are ignored.
• extractAllMatches(subject, mindex): Return all the matches in a single XStringViews object.

One could also use another function which does not give the position * match_fwd <- vcountPattern(fwd, seq,max.mismatch=0, min.mismatch=0, with.indels=FALSE, fixed=FALSE, algorithm=“auto”)

8.1 Summarize all eukaryotes

  pr2_match_summary_primer_set <- pr2_match_final_filtered  %>% 
                      group_by(gene_region, primer_label, primer_set_id ) %>% 
                      summarize (pct_fwd = sum(!is.na(fwd_pos))/n()*100,
                                 pct_rev = sum(!is.na(rev_pos))/n()*100,
                                 pct_amplicons = sum(!is.na(ampli_size))/n()*100, 
                                 ampli_size_mean = mean(ampli_size, na.rm=TRUE),
                                 ampli_size_sd = sd(ampli_size, na.rm=TRUE),
                                 ampli_size_max = max(ampli_size, na.rm=TRUE),
                                 ampli_size_min = min(ampli_size, na.rm=TRUE),
                                 n_seq = n())  
 
  write_tsv(pr2_match_summary_primer_set, str_c("pr2_match_summary_primer_set_",gene_selected,".tsv"), na="")
  
# Long form for Number of sequences
  
  # This dataframe is used to re-order the bars correctly
  pct_category_order <- data.frame(pct_category = c("pct_amplicons","pct_fwd","pct_rev"), pct_category_order = c(3, 1, 2))
  
  pr2_match_summary_primer_set_long <- pr2_match_summary_primer_set %>% 
    tidyr::gather("pct_category", "pct_seq", pct_fwd:pct_amplicons) %>% 
    left_join(pct_category_order)

8.2 Summarize per supergroup

  pr2_match_summary_primer_set_sg <- pr2_match_final_filtered  %>% 
                      group_by(supergroup, gene_region, primer_label, primer_set_id) %>% 
                      summarize (pct_fwd = sum(!is.na(fwd_pos))/n()*100,
                                 pct_rev = sum(!is.na(rev_pos))/n()*100,
                                 pct_amplicons = sum(!is.na(ampli_size))/n()*100, 
                                 ampli_size_mean = case_when (pct_amplicons>0 ~ mean(ampli_size, na.rm=TRUE)),
                                 ampli_size_sd = case_when (pct_amplicons>0 ~ sd(ampli_size, na.rm=TRUE)),
                                 ampli_size_max = case_when (pct_amplicons>0 ~ max(ampli_size, na.rm=TRUE)),
                                 ampli_size_min = case_when (pct_amplicons>0 ~ min(ampli_size, na.rm=TRUE)),
                                 n_seq = n())  %>% 
                      ungroup()
 

  write_tsv( pr2_match_summary_primer_set_sg, str_c("pr2_match_summary_primer_set_",gene_selected,"_per_sg.tsv"), na="")

8.3 Summarize per class

  pr2_match_summary_primer_set_class <- pr2_match_final_filtered %>% 
                      group_by(supergroup, division, class, gene_region, primer_label, primer_set_id) %>% 
                      summarize (pct_fwd = sum(!is.na(fwd_pos))/n()*100,
                                 pct_rev = sum(!is.na(rev_pos))/n()*100,
                                 pct_amplicons = sum(!is.na(ampli_size))/n()*100, 
                                 ampli_size_mean = case_when (pct_amplicons>0 ~ mean(ampli_size, na.rm=TRUE)),
                                 ampli_size_sd = case_when (pct_amplicons>0 ~ sd(ampli_size, na.rm=TRUE)),
                                 ampli_size_max = case_when (pct_amplicons>0 ~ max(ampli_size, na.rm=TRUE)),
                                 ampli_size_min = case_when (pct_amplicons>0 ~ min(ampli_size, na.rm=TRUE)),
                                 n_seq = n())  
 
  write_tsv(pr2_match_summary_primer_set_class, str_c("pr2_match_summary_primer_set_",gene_selected,"_per_class.tsv"), na="")  

9.1 Average size

  ggplot(pr2_match_final, aes(x= primer_label, y=ampli_size, group=primer_set_id)) + 
    geom_boxplot(outlier.alpha = 0.3) +
    theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    xlab("Primer set")

9.2 Plot an example of amplicon distribution (Fig. 3)

fig3A <- list()
fig3B <- list()


for (one_primer_set in c(4, 29)){
  
  if (one_primer_set == 4 ) {
    xmin=570
    xmax = 610
    xmax2 = 2000
  } else {
    xmin=130
    xmax = 190
    xmax2 = 1000
    
  }
  
  pr2_match_final_one <- pr2_match_final %>% 
    filter(primer_set_id==one_primer_set) %>% 
    filter( !(supergroup %in% c("Apusozoa", "Eukaryota_X", "Protoalveoalata") ))
    
  primer_label <- str_replace_all(pr2_match_final_one$primer_label[1], "_", " ")
  
  g <- ggplot(pr2_match_final_one, aes(x= ampli_size)) + 
    geom_density(fill="blue", alpha=0.9) +
    xlab("Amplicon size") +
    ggtitle(str_c("Primer set - ", primer_label)) +
    xlim(xmin,xmax)
  
  print(g)
  
  fig3A[[one_primer_set]] <- ggplot(pr2_match_final_one, aes(x= ampli_size, fill=supergroup)) + 
    geom_density(alpha=0.9)+
    theme_bw() +
    theme(legend.position = "none") +
    scale_fill_viridis_d(option = "inferno") +
    xlab("Amplicon size (bp)") + ylab("Density") +
    ggtitle(str_c("Primer set - ", primer_label)) +
    xlim(xmin,xmax) 
  
  print(fig3A[[one_primer_set]])

  fig3B[[one_primer_set]] <- ggplot(pr2_match_final_one, aes(x= supergroup, y=ampli_size)) + 
    geom_boxplot(outlier.alpha = 0.3)+
    theme_bw() +
    coord_flip()+
    # theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    xlab("Supergroup") + ylab("Amplicon size (bp)") +
    ylim(0,xmax2)
    
  print(fig3B[[one_primer_set]])
  
  g <- ggplot(pr2_match_final_one, aes(x= supergroup, y=ampli_size)) + 
    geom_violin() +
    coord_flip()+
    # theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    xlab("Supergroup") + ylab("Amplicon size (bp)") 
  print(g)
    
}

9.3 Plot an example of amplicon size distribution

fig3C <- list()

for (one_primer_set in c(4, 29)) {
  
  pr2_match_summary_filtered <- pr2_match_summary_primer_set_sg %>% 
    filter((n_seq > 20) &
           (primer_set_id == one_primer_set) & 
           !(supergroup %in% c("Apusozoa", "Eukaryota_X") ))
  

  
 fig3C[[one_primer_set]] <- ggplot(pr2_match_summary_filtered) + 
    geom_col(data=pr2_match_summary_filtered, 
             aes(x=str_c(supergroup, " - n = ", n_seq), 
                 y=pct_amplicons, 
                 fill=supergroup),  
                 position="dodge") +
    theme_bw()+
    theme(legend.position = "none") +
    scale_fill_viridis_d(option = "inferno") +
    coord_flip() + 
    ylab("% of sequences amplified") + xlab("") # +
    # ggtitle (str_c("Set - ", one_primer_set, " - % amplified per Supergroup") ) 

  print(fig3C[[one_primer_set]])  
  
  g <- ggplot(filter(pr2_match_summary_filtered, !is.nan(ampli_size_mean))) + 
    geom_point(aes(x=str_c(supergroup, " - n = ", n_seq), 
                   y=ampli_size_mean), 
               colour="black") +
    coord_flip()  +
    geom_errorbar(aes(x=str_c(supergroup, " - n = ", n_seq), 
                      ymax=ampli_size_mean + ampli_size_sd, 
                      ymin=ampli_size_mean - ampli_size_sd)) +
    xlab("Supergroup") + ylab("Amplicon size (bp)")   # +
    # scale_y_continuous(breaks = (1:8)*200, limits = c(0,1500)) +
    # ggtitle (str_c("Set - ", one_primer_set, 
    #                " - Amplicon sizes") )  +
    # geom_hline(yintercept = c(450,550) , linetype= 2)

  print(g)


  

}

9.4 Fig. 3

# g.fig3A <- ggplotGrob(fig3A) # convert to gtable
# g.fig3B <- ggplotGrob(fig3B) # convert to gtable
# 
# fig3A.widths <- g.fig3A$widths[1:3] # extract the first three widths, 
#                                   # corresponding to left margin, y lab, and y axis
# fig3B.widths <- g.fig3B$widths[1:3] # same for mpg plot
# max.widths <- grid::unit.pmax(fig3A.widths, fig3B.widths) # calculate maximum widths
# g.fig3A$widths[1:3] <- max.widths # assign max. widths to iris gtable
# g.fig3B$widths[1:3] <- max.widths # assign max widths to mpg gtable

cowplot::plot_grid(fig3A[[4]], fig3B[[4]], fig3C[[4]],fig3A[[29]], fig3B[[29]], fig3C[[29]], 
                   labels = c("A","","", "B","",""), ncol = 3, nrow = 2, align = "h")

10.1 All Eukaryotes

Comments

• Percent of sequences amplified
  • Logically, the % of seq amplified is always < min(% of sequences matching forwar, % of sequences matching reverse)
  • In general it is the reverse primer that causes problems.
  • Some primer sets do not amplify any sequence (11, 19, 20, 21)
  • For V9, primer set 30 reverse primer is in the ITS region which is not present in the PR2 sequences, so no amplification.
• Amplicon sizes
  • Only 8 V4 primer sets are suitable for Illumina 2x250 (max amplicon size = 450 bp)
  • This is if we consider the mean… If we consider the variation around the mean then, only 3 suitable for 2x250
  • Five more V4 primer sets are suitable for Illumina 3x250 (max amplicon size = 550 bp)

for (one_region in gene_regions) {
  
  pr2_match_summary_primer_set_region_long <- filter(pr2_match_summary_primer_set_long, gene_region== one_region)  
  pr2_match_summary_primer_set_region <- filter(pr2_match_summary_primer_set, gene_region== one_region)
  pr2_match_region <- filter(pr2_match_final, gene_region == one_region)
  
  g <- ggplot(pr2_match_summary_primer_set_region_long) + 
    geom_col(aes(x=reorder(primer_label, primer_set_id), y=pct_seq, 
                 fill=fct_reorder(pct_category, pct_category_order)), width=.7, position = "dodge") +
    theme_bw() +
    xlab("Primer set") + ylab("% of sequences amplified") + 
    scale_fill_manual(name = "% amplified",  values = c("pct_amplicons" = "black", "pct_fwd" = "blue","pct_rev" = "red"), 
                      labels=c( "Primer fwd", "Primer rev", "Amplicons")) +
    ggtitle (str_c(one_region, " - Percentage of sequences recovered"))+
    theme(axis.text.x = element_text(angle = 45, hjust = 1))
  
  print(g)


   g <- ggplot(filter(pr2_match_summary_primer_set_region, !is.nan(ampli_size_mean))) + 
    geom_point(aes(x=reorder(primer_label, primer_set_id), y=ampli_size_mean), colour="black") +
    geom_errorbar(aes(x=reorder(primer_label, primer_set_id), 
                      ymax=ampli_size_mean + ampli_size_sd, 
                      ymin=ampli_size_mean - ampli_size_sd)) +
    theme_bw() +
    xlab("Primer set") + ylab("Amplicon size (bp)") +
    # scale_y_continuous(breaks = (1:8)*200, limits = c(0,1500)) +
    ggtitle (str_c(one_region, " - Amplicon size - Lines correspond to limits for Illumina 2x250 and 2x300 respectively") ) +
    geom_hline(yintercept = c(450,550), linetype = c(2,3) )+
    theme(axis.text.x = element_text(angle = 45, hjust = 1))
   
  print(g)
  
     g <- ggplot(pr2_match_region) + 
    geom_boxplot(aes(x=reorder(primer_label, primer_set_id), y=ampli_size), colour="black", outlier.alpha = 0.3) +
    theme_bw() +
    xlab("Primer set") + ylab("Amplicon size (bp)") +
    # scale_y_continuous(breaks = (1:8)*200, limits = c(0,1500)) +
    ggtitle (str_c(one_region, " - Amplicon size - Lines correspond to limits for Illumina 2x250 and 2x300 respectively") ) +
    geom_hline(yintercept = c(450,550), linetype = c(2,3) )+
    theme(axis.text.x = element_text(angle = 45, hjust = 1))
   
  print(g)
    }

10.2 By supergroup

Comments

• Excavata have a very different patterns from the rest of the group. They are not amplified by quite a few primer sets. They have also bigger amplicons
• Some groups exhibit higher variability in amplicon size (e.g Chlorophyta)

for (one_region in gene_regions) {
  
  pr2_match_summary_primer_set_sg_region <- filter(pr2_match_summary_primer_set_sg, 
                                                   (gene_region == one_region) &(n_seq > 20)) 
  
  g <- ggplot(pr2_match_summary_primer_set_sg_region) + 
    geom_col(aes(x=supergroup, 
                 y=pct_amplicons),  
             fill="grey", 
             position="dodge") +
    theme_bw() +
    coord_flip() + 
    ylab("% of sequences amplified") + 
    xlab("Supergroup")  +
    ggtitle (str_c(one_region, " - % amplified per Supergroup") ) +
    facet_wrap(~ primer_label, scales ="fixed")

  print(g)  
  
  g <- ggplot(filter(pr2_match_summary_primer_set_sg_region, !is.nan(ampli_size_mean))) + 
    geom_point(aes(x=supergroup, 
                   y=ampli_size_mean), 
               colour="black") +
    theme_bw() +
    coord_flip()  +
    geom_errorbar(aes(x=supergroup, 
                      ymax=ampli_size_mean + ampli_size_sd, 
                      ymin=ampli_size_mean - ampli_size_sd)) +
    xlab("Supergroup") + 
    ylab("Amplicon size (bp)")  +
    ggtitle (str_c(one_region, " - Amplicon size - Lines correspond to limits for Illumina 2x250 and 2x300 respectively") ) +
    geom_hline(yintercept = c(450,550) , linetype= 2) +
    facet_wrap(~ primer_label, scales = "fixed")

  print(g)



}

10.3 By class for autotrophs

for (one_region in gene_regions) {

  pr2_match_summary_filtered <- filter(pr2_match_summary_primer_set_class,
                                      (n_seq > 20) &
                                      (division %in% c("Haptophyta", "Dinoflagellata", 
                                                       "Chlorophyta", "Ochrophyta", "Cryptophyta")) &
                                      (gene_region == one_region))
  
  if(nrow(pr2_match_summary_filtered) > 0 ) {
  
  g <- ggplot(pr2_match_summary_filtered) + 
    geom_col(data=pr2_match_summary_filtered, 
             aes(x=str_c(str_trunc(division, 3, ellipsis = ""), "-", class), 
                 y=pct_amplicons),  
             fill="grey", position="dodge") +
    theme_bw() +
    coord_flip() + 
    ylab("% of sequences amplified") + xlab("Class")  +
    ggtitle (str_c("Set -", one_primer_set, " - % amplified per Class") )  +
    facet_wrap(~ primer_label, scales ="fixed")

  print(g)
  
g <- ggplot(filter(pr2_match_summary_filtered, !is.nan(ampli_size_mean))) + 
    geom_point(aes(x=str_c(str_trunc(division, 3, ellipsis = ""), "-", class), 
                   y=ampli_size_mean), 
               colour="black") +
    theme_bw() +
    coord_flip()  +
    geom_errorbar(aes(x=str_c(str_trunc(division, 3, ellipsis = ""), "-", class), 
                      ymax=ampli_size_mean + ampli_size_sd, 
                      ymin=ampli_size_mean - ampli_size_sd)) +
    xlab("Class") + ylab("Amplicon size (bp)")  +
    # scale_y_continuous(breaks = (1:8)*200, limits = c(0,1500)) +
    ggtitle (str_c("Set -", one_primer_set, 
                   " - Amplicon size - Lines correspond to limits for Illumina 2x250 and 2x300 respectively") ) +
    geom_hline(yintercept = c(450,550) , linetype= 2) +
    facet_wrap(~ primer_label, scales ="fixed")

  print(g)
  }
}

11.1 Specific et sets for Opisthokonta

• 35 UnNonMet
• 16 Piredda
• 17 Comeau

for (one_primer_set in c(16, 17, 35)) {
  
  pr2_match_summary_filtered <- filter(pr2_match_summary_primer_set_class,
                                      (n_seq > 20) &
                                      (supergroup %in% c("Opisthokonta")) &
                                      (primer_set_id == one_primer_set))
  
  g <- ggplot(pr2_match_summary_filtered) + 
    geom_col(data=pr2_match_summary_filtered, 
             aes(x=str_c(division, "-", class, " - n= ", n_seq), 
                 y=pct_amplicons),  
             fill="grey", position="dodge") +
    theme_bw() +
    coord_flip() + 
    ylab("% of sequences amplified") + xlab("Class")  +
    ggtitle (str_c("Set -", one_primer_set, " - % amplified per Class") ) 

  print(g)  
  
  g <- ggplot(filter(pr2_match_summary_filtered, !is.nan(ampli_size_mean))) + 
    geom_point(aes(x=str_c(division, "-", class, " - n= ", n_seq), 
                   y=ampli_size_mean), 
               colour="black") +
    coord_flip()  +
    geom_errorbar(aes(x=str_c(division, "-", class, " - n= ", n_seq), 
                      ymax=ampli_size_mean + ampli_size_sd, 
                      ymin=ampli_size_mean - ampli_size_sd)) +
    xlab("Class") + ylab("Amplicon size (bp)")  +
    # scale_y_continuous(breaks = (1:8)*200, limits = c(0,1500)) +
    ggtitle (str_c("Set -", one_primer_set, 
                   " - Amplicon size - Lines correspond to limits for Illumina 2x250 and 2x300 respectively") ) +
    geom_hline(yintercept = c(450,550) , linetype= 2)

  print(g)

}